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ABSTRACT 

We develop a technique to investigate the possibility that some of the recently discovered ultra- 
faint dwarf satellites of the Milky Way might be cusp caustics rather than gravitationally self-bound 
systems. Such cusps can form when a stream of stars folds, creating a region where the projected 2-D 
surface density is enhanced. In this work, we construct a Poisson maximum likelihood test to compare 
the cusp and exponential models of any substructure on an equal footing. We apply the test to the 
Hercules dwarf (d ~ 113 kpc, My ~ —6.2, e ~ 0.67). The flattened exponential model is strongly 
favored over the cusp model in the case of Hercules, ruling out at high confidence that Hercules is a 
cusp catastrophe. This test can be applied to any of the Milky Way dwarfs, and more generally to 
the entire stellar halo population, to search for the cusp catastrophes that might be expected in an 
accreted stellar halo. 

Subject headings: Galaxy — halo; galaxies — dwarfs 



1. INTRODUCTION 

In the last five years, more than a dozen ultra-faint 
Milky Way satellites have been discovered with one mil- 
lionth of the Milky Way's luminosity or less. These ob- 
jects have stars with [Fe/ Hl at least as low as -3.0, [Fe/H] 
spreads of up to 1 dex (iKirbv et al.l l2008t iFrebel et~aLl 
12010b ISimon et al.l I2010D. and velocity dispersions of 
~ 3 - 7 km sec -1 (|Simon fc G"ehal 120071: IMartin et all 
I2007D . and scale sizes of 30 - 200 pc (jMartin et al.ll2008t) . 
A common explanation for these observed properties thus 
far is that we are observing dwarf galaxy member stars 
moving on bound orbits within highly dark matter dom- 
inated potential wells. Within this scenario, the veloc- 
ity dispersions of these objects have been interpreted to 
possibly imply mass-to- light ratios of 10 2 ~ 4 within their 
central 300 pc (jStrigari et al.ll2008h and that galaxies can 
have total luminosities of only ~ IOOOLq - far less than 
the luminosities of some individual stars. 

These remarkable conclusions should require a high 
burden of proof: Are there any other plausible expla- 
nations for the ultra-faint Milky Way satellites as a pop- 
ulation? One possibility is that the ultra-faints could 
have intrinsic near-zero velocity dispersions (and thus no 
dark matter), but that binary stars have inflated their ob- 
served dispersions. While this scenario could be conceiv- 
able for a small number of objects, it cannot explain away 
the entire ultra-faint galaxy population (ISimon et al.l 
l20lfi iMartlnez et al.ll20iol iMcConnachie fc C6tell2010h . 
Moreover, the low amount of and large spread in their 
iron abundances show that these objects cannot have 
formed in the same way as globular clusters, ruling 
out (disrupting) globular cluster scenarios. The possi- 
ble inflated ellipticities of ultr a-faint dwarfs relat ive to 
their more luminous neighbors (jMartin et al.|[2008h com- 
bined with claims of irregular isophotes tempts a con- 
clusion that these are tidal remnants. However, claims 
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of isophote irregularity have been inconc lusive for the 
major ity of ultra- faint satellites. Recently. IMartin fe; Jin! 
(2010) have found that the shape of the Hercules dwarf 
might be explained by tidal disruption. 

A fourth possibility is that (a subset of) the kinemat- 
ically cold, spatial overdensities that have been identi- 
fied as ultra-faint galaxies are instead cusp caustics of 
cold stellar debris . Cold substru c tures a re both observed 
(iGrillmairl 120061 : lYannv et al.l 120031 : iBelokurov et~aTl 
120061 e.g?) and predicted to be abundant in the Milky 
Way's halo. Shell-like structures can form during the 
disruption of a low mass satellite, as the satellite inter- 
acts with a massive galaxy, and have been seen in both 
simulations (iHernquist fc Quinnlll988l) and observations 
(jMalin fc Carterlll983l ) of early- type galaxies. After an 
interaction creates such a sheet of stars, the sheet can 
evolve and fold in ways that produce visible fold caustics 
("catastrophes"), cusp caustics, and higher-order caus- 
tics in the observ ed two-dimensional distribution of stars 
()TremaindH999T ). A cusp caustic in a two-dimensional 
observable space appears as a highly elliptical centroid 
of stars where two fold caustics meet, creating an asym- 
metrical flare-like feature. In this scenario, some of the 
ultra-faint dwarfs may actually be cusps and some of the 
cold streams could be folds. 

In this paper, we develop a technique for testing the 
hypothesis that ultra-faint dwarf galaxies are the 2-D 
projections of such folded stellar sheets. We apply this 
technique to the Hercules Milky Way satellite , which has 
d ~ 113 kpc, M v ~ -6.2 , and e - 0.67 (|Sand et al.l 
120091: IColeman et ail l2007t) . We choose this object be- 
cause it has the highest ellipticity of any of the Milky 
Way's dwarfs and is thus a good candidate for the cusp 
model. 

2. OBSERVATIONAL DATA 

IColeman et al.l (|2007l hereafter C07) obtained Large 
Binocular Telescope imaging of the Hercules dwarf 
galaxy in the Gunn r band (25 minutes), V band (20 
minutes), and B band (30 minutes) to derive its struc- 
tural parameters. The Large Binocular Camera spans a 
23' x 23' field of view. C07 applied color-color-magnitude 
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filtering in colors c% and C2 (combinations of r, V, and 
B) to their point source catalog to select objects most 
consistent with belonging to Hercules rather than to the 
field. In this paper we utilize the ~10,000 stars from 
the C07 CMD selected subset of their LBT/LBC Her- 
cules catalog. This catalog is identical to what they used 
to derive the structure of Hercules. Please see C07 for 
further details on this data set. 

3. MODEL COMPARISON: CUSP VS EXPONENTIAL 

We aim to determine whether a cusp model is a better 
explanation of the spatial distribution of the elongated 
Hercules dwarf than a flattened exponential model. To 
perform this model comparison, we specify each model 
explicitly, compute the Poisson likelihood of the data 
given the models, and use the maximum likelihood of 
each model to determine which best fits the data. 

3.1. Constructing the Models 

Both of the exponential and cusp models tested here 
have six free parameters. The cusp model's free parame- 
ters are: central RA and Dec, rotation angle and scale, a 
smoothing scale, and a constant background level. The 
exponential model's free parameters are: central RA and 
Dec, a rotation angle, ellipticity, half-light radius, and 
a constant background level. For each model we step 
through the six parameters (0) and evaluate the likeli- 
hood of the data given each model. The full likelihood, 
which is a product of the individual star likelihoods, is: 

^(®) = n ^0) s ^i ) w 

i—l * ' 
^ Nrandom 

Z(&) = — ^random,j\®) (2) 

^random . 

3 = 1 

where is the list of parameters for the cusp or exponen- 
tial model, £(•) is the surface density of the model given 
the parameters, and Z(®) is a normalization constant, 
estimated using a spatially randomized stellar catalog 
over the observed sky region. N sta r is the total number 
of stars, and x, is the two-dimensional position (on the 
sky) of star i. N ran( i om is the total number of random 
points, set to equal N star , and x ran( iom.j is the position 
of random point j . 

3.1.1. The Cusp Model 

The cusp surface density m odel is computed following 
the work of iTremaind (|1999t h First, a background- free 
S cusp is computed on a grid in the natural coordinate 
system £: 

W6,6)=£[1 + t^|^] 1/2 , (3) 

roots ^ 2 y > 

where 

a = -y 3 - by- (4) 

In three dimensions, coordinates £i and £2 make up the 
2-D observable space, while coordinate y is the third, 
and hidden, dimension. As a third degree polynomial, 
equation (4) has either one or three solutions, depending 



on the position in £. Equation (3) is to be summed over 
all the possible roots of variable y, obtained by solving 
equation (4). In regions where there is no cusp, there 
will only be one solution to equation (4), and hence no 
summation in the calculation of the surface density is 
necessary. On the other hand, in regions of £ where there 
is a cusp, the variable y will be triple valued, and equa- 
tion (3) must be summed over all three roots. Because 
this function has a large dynamic range, we compute it 
on a fine grid of 1000 x 1000 pixels. Even with such 
fine pixels, however, the value of the cusp function can 
vary significantly within a pixel in regions of the image 
where the cusp is present. At such pixels, we evaluate 
the function (equation 3) at multiple random locations 
within the pixel. This adaptive sampling sets the num- 
ber of samples per pixel to a constant multiple (10) of 
the value of the function. This essentially ensures that 
we are sampling more in the most interesting regions of 
the image, where the cusp is present. T, cusp has no free 
parameters in the natural coordinate system. 

After the cusp image has been computed on a grid, 
it is smoothed with a symmetric two-dimensional Gaus- 
sian kernel with variance a 2 . We have computed seven 
cusp models, each with a different smoothing scale, with 
values of a = [1, 2, 4, 8, 16, 32, 64] pixels. The two dimen- 
sional positions Xi (RA,Dec) of each star i on the sky 
is converted to a position £i in the natural coordinate 
system of the cusp model by a shift, scale and rotation: 

fc = R.[x,-X], (5) 

where R is a 2 x 2 matrix encoding the two free pa- 
rameters of scale and rotation, and X represents the two 
free parameters that shift the position of the center of 
the model on the sky. Finally the surface density £ com- 
puted at point £i is the sum of the background free T, cusp , 
and a constant background level £b ff : 

S(0 = £c« sp (0 + S fc5 . (6) 

3.1.2. The Exponential Model 

We follow the same procedure for the exponential 
model as for the cusp model. The exponential model, 
in its natural units, is: 

E«xp(0 = exp(-\£\), (7) 

and the model is not smoothed at all. The 2x2 matrix 
R used in this case is a symmetric tensor of scale and 
shear with three free parameters, used to calculate the 
position angle (8), ellipticity (e), and scale radius (r^) of 
the model. 

All models have been defined on grids of 1000 x 1000 
pixels, on the interval £=[-5,5], and have six degrees of 
freedom. Figure 1 shows the surface brightness grids 
for cusp models with a = [1,2,4,8] pixels and for the 
exponential model, each in their natural units. 

3.2. MCMC 

We use a Markov Chain Monte Carlo (MCMC) method 
to calculate the maximum likelihood of each model, and 
its corresponding parameters. MCMC simulates the like- 
lihood for a set of parameters by sampling from the pos- 
terior distribution through a series of random draws. 
We specifically use the Metropolis-Hastings Algorithm, 
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Maximum Likelihood Parameters 
Parameter Measured Uncertainty | lnjgf 



Best Exponential Profile 



R.A. (h m s) 16:31:03.9 8" 
DEC (d m s) +12:47:10 7" 
rh (arcmin) 5.3 0.4 
e 0.65 0.03 
9 (degrees) -75.6 2.1 
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Best Cusp Profile 


R.A. (h m s) 16:31:30 17" 
DEC (d m s) +12:46:15 4" 
scale a 37.6 2.1 
6 (degrees) -78.1 3.1 
cr (pixels) 2 b 
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which generates a random walk and steps through each 
parameter in such a way that more probable values of pa- 
rameter space are more often stepped to. For both the 
cusp and exponential models we choose a separable prior 
that is flat in center position (X) , flat in the components 
of the R matrix, and flat in the background level 
with a lower limit of zero. There were no derivatives in 
our prior at the peaks in the likelihood, so the choice of 
prior does not affect any of our results. The acceptance 
ratios for all parameters are ~ 0.5. Though less computa- 
tionally efficient, this method is similar to the ma ximum 
likelihood technique used by (jMartin et a.1.1 [2008h to de- 
termine the structural parameters of dwarf galaxies, and 
is better at escaping local minima. 

3.3. Results 

We find that a flattened exponential model is a far bet- 
ter fit to the Hercules data than any of the seven cusp 
models tested in this work. Of the seven cusp models 
tested, the cusp with a smoothing scale of a = 2 pixels 
(0.75 arcmin) was found to be the most probable. The 
difference between the maximum likelihood of the expo- 
nential and best fit cusp model is 

Nominally this rules out the best cusp model in favor of 
the exponential at a confidence level of > 99.99%. 

Table 1 lists the parameters at maximum likelihood for 
the exponential and best cusp models. The structural pa- 
rameters for the exponential model yield a position angle 
for Hercules at —76°, and an ellipticity of 0.65 (defined 
as e = I — b/a), both i n goo d agreement w i th the find- 
ings of ICo eman et all (|2007j ). iMartin et all (|2008t ). and 
ISand et "al (|2009f ). The parameters that maximize the 
likelihood of the cusp model are such that the model is 
slightly smoothed, and highly scaled. This is not sur- 
prising since the cusp model is a poor fit to the data, 
and at its optimal configuration it will be smoothed and 
scaled in such a way to resemble an flattened exponential- 
like distribution near the center. If Hercules were a 2-D 
projection of a cusp, the spatial distribution of its stars 
would show a highly elongated center of stars, along with 
two "tails" indicative of the two fold catastrophes. Fig- 
ure 2 shows the predicted spatial distribution of stars for 
the best fit exponential and cusp models. 

While the cusp model is strongly ruled out in favor of 
a flattened exponential, we make no claims that the ex- 
ponential profile is the best available model for Hercules. 

4. CONCLUSION & DISCUSSION 



It has been suggested that some of the ultra- faint dwarf 
satellites of the Milky Way may be the 2-D projection 
of cusp caustics, and not gravitationally bound galax- 
ies. The highly elliptical Hercules dwarf is a good candi- 
date for testing the cusp model, which predicts a highly 
elongated centroid of stars. We have established a hy- 
pothesis test in order to determine whether the Hercules 
overdensity is better modeled with a cusp or an exponen- 
tial, and rule out at very high confidence the possibility 
that Hercules is a cusp catastrophe. While other surface 
brightness m odels have been fit to Her cules data, such 
as the King (|Kingl ll966) and Plummer (|Plummerlll911| ) 
profiles, we have chosen to not test such models here, and 
so have not evaluated how they compare to the flattened 
exponential model. This work also has not evaluated if 
the properties of the Hercules galaxy are better explained 
by the exponential e quilibrium model th an by models of 
tidal disruption (e.g lMartin fc Jinll20I0[ ). Our intention 
in this work is not to advocate that an exponential pro- 
file is the best explanation of the observed ellipticity of 
Hercules, but rather to definitively remove from consid- 
eration the possibility that Hercules is a cusp caustic. 

The work outlined in this paper can be used to test 
any of the Milky Way dwarfs suspected of being cusps. 
The highly elliptical dwarfs Ursa Major I and II are in- 
teresting candidates for such a test, though both lack 
obvious evidence of the double spatial tail indicative 
of the fold caustics that meet at the cusp. Further- 
more, the cusp caustic predicts triple-valued velocities 
for stellar members of the cusp. Spectroscopic data of 
the ultra-faint dwarfs suspected of being cusps can thus 
be used to further rule out or confirm such hypotheses. 
For example, spectroscop ic data in the vicinity of the 
ultra- faint galaxy Segue I (jNiederste-Ostholt et al.ll20"oll 
ISimon et al.1 l2010l) show a lot of interesting structure, 
and could be used to test if Segue I is a cusp. 

The techniques introduced in this paper can be fur- 
ther extended to implement a search for cusp and fold 
caustics on the entire sky, using data from a large sur- 
vey like the Sloan Digital Sky Survey (SDSS), Gaia, 
or the Large Synoptic Survey Telescope (LSST). Stel- 
lar streams and dwarf galaxies, which can be thought 
of, respectively, as one and (blurry) zero dimensional 
objects in phase phase, have been found with multi- 
band two-dimensional angular maps in photometric sur- 
veys, such as SDSS. Such data sets are two-dimensional, 
though some distance information can be determined at 
low signal-to-noise using colors, and are well suited for 
finding low dimensional structures. Given the number 
of streams already known in the Milky Way's stellar 
halo, it is likely that higher dimensional structures in 
six-dimensional phase s pace exist, but have not yet been 
found in current data (|Tremainelll999l ). Future surveys 
that return full (or nearly full) phase-space info rmation, 
such as the Gaia mission (jPerryman et al.ll200fl ). will be 
much better suited than current surveys for finding sys- 
tems with higher dimensionality. Meanwhile, such struc- 
tures may only be easily identifiable in photometric data 
if their projection into 2-D observable space creates re- 
gions of enhanced density, like the cusp and fold caustics. 



We thank Niayeshi Afshordi for inspiring this paper, 
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Fig. 1. — From left to right: The surface brightness maps produced by cusp models smoothed with a Gaussian kernel of variance 
a = [1,2,4,8] pixels and by an exponential model, as described in Section 3.1. Each was calculated on a 1000 x 1000 pixel grid, and is 
plotted in its natural coordinate system, §. 
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Fig. 2. — Left Panel- The observed data set of Hercules stars from [Coleman et alJ l|2007fl . The elongated spatial distribution of this 
overdensity has raised questions as to whether it is a bound dwarf galaxy or a cusp fold. Middle Panel- The predicted spatial distribution 
of stars in the exponential model. Right Panel- The predicted spatial distribution of stars in the cusp model with the highest likelihood 
(smoothed with a = 2 pixels). The dashed red lines mark the best fit exponential and cusp models in middle and left panels, respectively. 
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